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1.  INTRODUCTION  AND  BACKGROUND 

Improved  3D  wave  speed  models  of  the  greater  Middle  East,  from  the  Turkish  Plateau  to  the 
eastern  edge  of  Tibet  (Figure  1),  enhances  our  ability  to  discriminate  between  natural  and  man¬ 
made  events,  locate  these  events,  identify  source  depths,  and  detennine  magnitudes  using 
wavefonn  based  methodologies.  Current  wave-speed  models  of  the  Middle  East  are  improved 
through  an  adjoint  tomography  method  [1]  applied  to  seismic  phase  and  amplitude  data.  An 
initial  step  towards  the  adjoint  tomography  model  requires  the  development  of  a  database  of 
relocated  and  well-characterized  sources  and  wavefonns.  The  adjoint  tomography  method  [1] 
uses  full  seismic  waveforms  as  a  measure  of  misfit  of  the  current  model  iteration.  Differences 
between  data  and  synthetics  are  used  to  create  adjoint  sources  and  generate  sensitivity  kernels 
required  to  update  the  wave-speed  model.  The  adjoint  tomography  model  was  built  on  the 
foundation  of  the  seismic  waveform  database  and  finite-frequency  wave  propagation,  and  will  be 
distributed  to  the  community 

Improvements  to  the  3D  wave  speed  model  for  the  Middle  East  using  the  adjoint  tomography 
method,  exploits  the  full  seismic  waveform,  including  amplitude  and  phase.  Due  to  these 
requirements  of  improving  an  existing  3D  model  and  using  seismic  waveform  data,  synthetic 
data  need  to  be  accurate  for  the  chosen  3D  wave  speed  model  and  the  sources,  e.g.  initial  event 
conditions  need  to  be  precise  when  comparing  to  real  sources  or  events. 

To  assess  the  precision  of  our  synthetic  sources,  events  were  re-inverted  in  ID  and  3D  wave- 
speed  models  and  agree  fairly  well  with  currently  available  solutions.  Use  of  a  3D  model  further 
improves  the  wavefonn  fit  between  data  and  synthetics,  reduces  the  overall  enor  of  the  CMT 
solution,  provides  a  means  to  include  shorter  period  arrivals,  and  improves  most  metrics  of 
wavefonn  fit. 

Following  source  re-interpretation,  the  3D  wave  speed  model  was  updated  using  adjoint 
tomographic  methods.  The  power  of  the  methodology  allows  for  incorporation  of  true,  physics- 
based,  wavefield  propagation  to  localize  errors  in  the  wave  speed  model.  Once  errors  are 
localized  and  identified  the  model  can  be  updated  and  another  iteration  is  possible.  Final 
iterations  for  the  Middle  East  3D  wave  speed  model  amplify  and  sharpen  wave-speed  anomalies 
while  increasing  the  predictive  capabilities  at  periods  of  20  seconds  and  longer  for  surface  waves 
and  14  seconds  and  longer  for  body  waves.  This  results  in  a  wave-speed  model  with  a  resolution 
of  about  60  km. 
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Figure  1.  Map  of  Re-Characterized  Events  in  the  Middle  East 


2.  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

The  adjoint  tomographic  inversion  is  a  two-step  procedure  where  a  source  re-interpretation  is 
first  perfonned  to  develop  a  more  precise  realization  of  the  events  within  the  region.  This  is 
followed  by  an  iterative  3D  wave  speed  inversion  driven  by  adjoint  methods. 
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2.1.  Event  Re-Characterization 


Within  the  Middle  East  approximately  200  events  were  re-characterized  using  the 
Moment  Tensor  (MT)  inversion  methodology  of  [2]  and  both  the  ID  PREM  model  [3]  and  a  3D 
wave-speed  model  (S2.9EA)  [4]  at  periods  of  25-125  seconds  and  of  60-125  seconds.  A  map  of 
all  available  stations  and  events  are  shown  in  Figure  1 .  Initially,  the  procedures  for  processing 
large  amounts  of  data,  comparing  the  data  to  synthetics,  and  reevaluating  the  source  parameters 
and  locations  needed  to  be  assessed  and  streamlined.  This  development  was  accomplished 
through  the  use  of  a  ID  wave-speed  model  and  long-period  modeling,  and  needed  to  be 
straightforward  and  nimble  enough  to  avoid  problematic  areas,  as  the  initial  3D  event  re- 
evaluation  requires  more  demanding  computational  efforts.  Problematic  events  and  stations  were 
subsequently  screened  out  within  this  procedure  before  building  the  final  waveform  dataset. 
Problematic  stations  and  waveforms  with  dropouts  and  poorly  characterized  amplitude  responses 
can  negatively  influence  a  MT  inversion,  and  these  stations  were  removed  before  any  3D  MT 
inversion  or  adjoint  inversions  were  perfonned.  Along  with  the  development  of  a  semi- 
automated  procedure,  the  data  were  evaluated  and  bad  quality  data  were  omitted  from  the  entire 
source  and  wave  speed  improvement  procedure.  The  full  inversion  procedure  we  followed  is 
detailed  in  [5]  including  the  evaluation  function  and  efforts  to  assess  the  solution  stability. 


1997/05/13 

II.KURK.LHZ  ( 36.51, 70.68)  |u_K|EV1HZ 


0  200  400  600  800  1000  200  400  600  800  1000  1200  1400 


seconds  seconds 


2005/02/22 

II.KURK.LHZ  <30-76-5681  >  IU.KIEV.LHZ 


200  400  600  800  1000  1200  200  400  600  800  1000  1200  1400 

seconds  seconds 


Figure  2.  Comparison  Between  ID  and  3D  Wave  Speed  Models  at  Long  and  Short  Periods 


Figure  2  demonstrates  the  need  to  use  an  initial  3D  wave  speed  model  as  opposed  to  a  ID 
model  as  the  starting  point  for  the  adjoint  inversions.  While  the  ID  model  was  used  as  a  testbed 
to  develop  the  procedures  to  accomplish  the  MT  inversion,  the  synthetics  and  their  prediction  of 
the  data  at  long  and  short  periods  is  important.  At  long  periods,  60-125  s,  top  row  of  Figure  2, 
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the  ID  and  3D  (not  shown)  wave  speed  models  perfonn  equally  well  at  data  reproduction.  The 
primary  arrivals  at  this  period  are  the  fundamental  mode  surface  wave  and  very  long  period  body 
waves.  However,  using  the  same  model,  synthetics,  and  data  but  at  a  shorter  period,  middle  row 
of  Figure  2,  highlights  the  deficiencies  in  the  ID  wave  speed  model.  Again,  the  fundamental 
mode  of  the  surface  wave,  the  largest  amplitude  and  longest  period  arrival,  is  well  fit  by  the  ID 
wave  speed  model,  but  higher  mode  and  shorter  period  surface  waves  are  completely 
unaccounted  for  by  the  ID  model.  Event  2005/02/22,  right  columns  of  Figure  2,  demonstrates 
the  lack  of  predicting  power  for  the  shorter  period  surface  waves  highlighted  in  the  blue  boxes. 

In  contrast,  the  3D  wave  speed  model  of  the  same  events  does  much  better  at  matching  the 
dispersion  and  reproducing  the  amplitudes  of  the  surface  waves.  Body  waves  at  shorter  periods 
are  also  matched  much  better  by  the  3D  wave  speed  model.  The  left  columns  of  Figure  2,  for 
event  1997/05/13,  show  that  synthetics  based  on  the  3D  model  better  fit  the  body  wave  (blue 
boxes)  phases  and  amplitudes  than  do  synthetics  based  on  the  ID  model. 
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Figure  3.  Comparison  of  the  MT  Solutions  From  Initial  to  Final  Solutions 


As  stated  above,  the  use  of  a  3D  model  at  shorter  period  is  essential  for  the  adjoint 
inversion.  This  is  due  to  a  major  reduction  in  computational  cost  as  the  a  3D  initial  model  is 
about  7  adjoint  iterations  further  along  than  a  ID  initial  model.  Inversions  for  the  MT  also 
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benefit  greatly  from  the  vastly  improved  3D  model  because  1)  the  inversion  is  able  to  use  shorter 
period  data  and  2)  use  much  more  of  the  full  waveform  than  possible  with  a  ID  model.  Use  of 
better-predicted,  shorter  period  data  is  important  for  MT  inversions  especially  for  stations  near 
nodal  planes  of  the  MT.  Amplitudes  of  arrivals  near  the  nodal  planes  can  be  quite  small,  can 
change  rapidly,  and  have  the  largest  impact  (largest  derivatives)  on  the  MT  solution,  as  such  the 
shorter  period  data  allows  the  exploitation  of  these  smaller  amplitudes.  As  seen  in  Figure  3,  our 
solutions  are  quite  close  to  the  original  solutions  derived  from  long  period,  teleseismic  arrivals, 
but  use  of  the  shorter  period  data  and  nodal  stations  allows  for  a  large  variance  reduction 
improvement  over  the  original  MT  solutions.  Minor  changes  in  the  MT  solution,  Figure  3,  can 
have  large  impacts  on  nodal  stations  that  ultimately  drive  large  changes  in  variance  reduction. 
Figure  3  includes  three  different  MT  solutions,  which  demonstrate  our  ability  to  recover  the  MT 
solution  with  our  methodology. 

Full  seismic  waveform  data  was  used  to  accomplish  the  MT  inversion  with  time  window 
selection  from  the  FLEXWIN  tool  [6]  based  on  a  set  of  relaxed  constraints,  detailed  in  [5],  A 
variety  of  inversions  were  performed  to  assess  the  impact  of  station  distribution  and  the  level  of 
constraint  on  the  final  solution.  While  [5]  reports  solutions  for  6  independent  seismic  moment 
tensor  elements  (Mij,  zero-trace)  plus  depth,  solutions  using  weaker  or  stronger  constraints  are 
consistent  across  inversion  parameterization.  To  assess  the  errors  in  the  inversion,  we 
implemented  a  bootstrap  procedure  that  randomly  selects  the  different  seismic  stations  to  include 
in  the  inversion.  As  the  inversion  is  quite  fast,  we  are  able  to  easily  implement  the  bootstrap  at 
minimal  computational  cost  and  obtain  estimates  of  the  standard  error  of  the  solutions.  To 
evaluate  the  stability  and  consistency  of  each  event,  the  P-axis  of  the  MT  is  used  as  it  provides  a 
set  of  two  numbers,  i.e.  trend  and  plunge,  rather  than  six  numbers,  Mij,  to  compare.  To  note,  the 
T-axis  of  the  MT  shows  similar  behavior  to  the  P-axis  allowing  this  reduction. 

2.2.  Adjoint  Tomographic  Inversion 

Using  the  re-characterized  events  and  seismic  wavefonn  data  from  Section  2.1,  the  wave 
speed  model  was  improved  using  an  adjoint  tomographic  inversion  procedure.  The  inversion 
procedure  is  iterative  so  substantial  improvement  over  traditional  inversion  procedures  is 
possible  and  also  allows  continued  improvement  of  a  3D  wave  speed  model.  The  adjoint 
inversion  is  built  upon  the  ability  to  localize  errors  within  a  wave  speed  model  that  manifest  as 
mismatches  between  real  data  and  synthetic  seismograms  from  the  current  iteration  wave  speed 
model.  The  inversion  procedure  is  a  four  step  process:  1)  create  synthetic  seismograms  for  the 
current  wave  speed  model  [7],  2)  create  adjoint  sources  from  metrics  between  the  data  and 
synthetics  [1],  3)  use  adjoint  sources  and  the  back-propagated,  synthetic  wavefield  from  step  1  to 
localize  and  quantify  errors  in  the  model  [8],  4)  update  the  model  by  minimizing  the  errors  in  the 
model  and  continue  the  iterations.  Updates  to  the  Middle  East  model  were  focused  on  improving 
the  wave  speeds  and  the  adjoint  sources  were  constructed  to  emphasize  this. 

Synthetics  seismograms  are  generated  in  step  1)  using  the  SPECFEM3D  global  package  [7]. 
This  earthquake  simulation  tool  is  able  to  generate  synthetic  seismograms  from  3D  wave  speed, 
attenuation,  and  anisotropic  models  including  effects  of  gravitation,  rotation,  topography,  and 
bathymetry.  Included  within  the  package  is  the  ability  to  save  the  final  set  of  time-steps 
necessary  to  accomplish  the  back-propagation  or  time  reversal  in  step  3). 
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Figure  4.  Single  Station-Event  Frechet  Kernel  Computed  in  a  3D  Wave  Speed  Model 


To  create  the  adjoint  sources  in  step  2),  automatic  measurements  are  made  between  the  data 
and  synthetics  from  SPECFEM3D  using  the  FLEXWIN  tool  [6],  creating  data  windows  and 
associated  metrics.  Windows,  and  their  adjoint  sources,  are  not  constrained  to  any  specific 
seismic  phase  or  time  window  and  only  use  well  matched  data-synthetic  pairs.  In  step  3)  adjoint 
sources  (step  2),  a  time  series,  along  with  the  reconstructed,  full  synthetic  wavefield  (step  1)  are 
propagated  in  reverse  time  through  the  current  iteration  wave-speed  model  to  identify  locations 
in  the  model  requiring  improvement.  Figure  4  shows  an  example  Frechet  kernel  derived  from 
mismatches  between  data  and  synthetic  seismograms  from  one  event  to  one  station,  using  two 
horizontal  and  the  single  vertical  component.  For  reference,  Frechet  kernels  derived  from  ID 
models  resemble  simple  banana-doughnut  structures  for  longer  period  body  waves.  The  kernels 
displayed  in  Figure  4  are  more  complex  due  to  the  presence  of  3D  structures  and  the  use  of  the 
full  wavefield,  including  body  and  surface  waves.  Interactions  between  the  adjoint  sources  and 
the  time-reversed  wavefield,  integrated  over  time,  generate  these  kernels  specific  to  each 
measurement,  e.g.  a  simple  difference  between  the  data  and  synthetics.  Such  kernels  can  be 
created  using  individual,  isolated  measurements  (Figure  4)  or  for  each  event  with  a  large  number 
of  back-propagated  adjoint  sources  (Figure  5). 

The  kernels  derived  using  full-waveforms  and  a  3D  wave  speed  model  still  resemble  the 
traditional  banana-doughnut  shape,  but  have  been  modified  by  the  interactions  of  multiple 
Frechet  kernels  from  individual  arrivals,  and  much  more  severely  because  of  the  use  of  a  3D 
wave  speed  model.  Appearance  of  a  number  of  “caustics”  in  the  single  event-single  station 
kernel  in  Figure  4  are  caused  by  either  sharp  changes  in  the  current  wave  speed  model  or  more 
likely  structures  existing  in  the  Earth  but  unaccounted  for  in  the  3D  wave  speed  model.  Large 
values  in  the  kernel  such  as  the  Arabian  Sea  south  of  the  event  and  the  southern  extent  of  the 
Caspian  Sea  are  a  result  of  large  amplitude  arrivals  in  the  data  that  are  missing  in  the  synthetics. 
This  mismatch  is  emphasized  in  the  construction  of  the  adjoint  sources  that  are  then  back- 
propagated  with  the  original  wavefield  to  identify  these  two  locations.  The  Caspian  Sea  kernel 
excursion  is  most  likely  due  to  thickened  sediments,  >  12km,  in  the  southern  Sea,  while  the 
southern  Arabian  Sea  anomaly  is  a  reflection  from  the  ocean-continent  boundary.  Moreover, 
larger  scale  anomalies  spread  or  smeared  over  the  Frechet  kernel  normally  indicate  a  overall 
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wave  speed  mismatch  between  the  true  model,  the  real  Earth,  and  the  current  3D  wave  speed 
model. 

Event  kernels,  similar  to  Figure  4  but  with  more  stations,  are  then  summed  volumetrically  to 
produce  Frechet  derivates,  Figure  5,  that  are  used  to  update  the  wave-speed  model  to  complete 
step  4  and  a  single  iteration.  The  process  is  then  restarted,  or  iterated  on,  to  continually  improve 
the  3D  model.  As  a  consistency  check,  a  small  subset  of  events,  ~5%,  are  not  included  within 
the  inversion  process,  but  synthetics  are  calculated  and  metrics  for  the  evaluation  windows  are 
derived.  Figure  5  shows  5  event  kernels  summed  together.  The  individual  events’  kernels  are 
identifiable  in  the  right  side  of  Figure  5  as  large  values  in  the  Frechet  derivative.  A  good 
example  of  the  large  values  near  the  source  is  the  event  to  the  east  of  the  southern  Caspian  Sea. 
The  event  creates  an  inverted  “U”  shape  on  the  surface  surrounding  the  source. 


Figure  5.  Multiple  Station-Event  Frechet  Kernel  Computed  in  a  3D  Wave  Speed  Model 


3.  RESULTS  AND  DISCUSSION 
3.1.  Event  Re-characterization 

The  comparison  between  Global  CMT  solutions  and  the  MT  solutions  presented  here  show 
negligible  differences  in  source  mechanisms  as  measured  from  the  P-axis.  This  is  true  for  all 
cases  (ID  vs  3D  model,  minimum  period  25s  vs  60s)  examined  and  when  varying  constraints  are 
applied  to  the  inversion  [5].  In  all  cases,  a  reduction  in  variance  between  data  and  synthetic  is 
seen  between  the  original  solution  and  our  calculated  solutions.  The  greatest  reduction  in 
variance  is  seen  in  the  3D  case  when  filtered  between  25- 125s.  A  majority  of  events  see  a 
variance  reduction  of  5-40%,  which  represents  a  significant  improvement  in  fitting  waveforms. 
The  largest  variance  reductions  are  a  result  of  a  better  amplitude  agreement  between  the  data  and 
synthetics. 

To  assess  the  robustness  of  the  MT  solutions  we  implemented  a  bootstrap  procedure  to 
detennine  the  solutions’  standard  errors.  We  use  the  orientation,  trend,  and  plunge  of  a  MT’s  P- 
axis  to  quantify  its  error.  For  all  of  our  events,  the  P-axis  distributions  are  similar  to  the  T-axis 
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distributions.  For  two  events  the  trend  of  the  P-axis  is  shown  in  Figure  6  where  results  from  the 
ID  wave  speed  model  are  shown  as  an  outlined  histogram  and  the  3D  wave  speed  model  is  the 
gray  histogram.  Also  included  in  Figure  6,  are  the  lower  hemisphere  projections  of  the  P-axes 
for  all  bootstrap  runs  for  the  selected  events.  A  tighter  cluster  of  points  indicates  a  smaller  error 
in  trend  and  plunge  of  the  P-axis  and  of  the  MT  solution  overall.  For  all  events,  the  P-axis  has  a 
larger  standard  error  in  trend  than  in  plunge.  This  can  be  seen  as  the  cluster  of  points  is  normally 
smeared  along  circles  tangential  to  the  bounding  circle,  see  Event  1997/05/13  ID.  P-axes  were 
corrected  for  0/360  crossing  and  plunges  near  horizontal  using  a  simple  analysis  to  find  the 
smallest  standard  error  by  forcing  the  P-axes  to  be  within  a  similar  reference  frame. 


(a)  1997/05/13 


(b)  2005/02/22 
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Figure  6.  Results  of  the  Bootstrap  Procedure  to  Determine  Standard  Errors  for  Two 

Events 


The  largest  contribution  to  the  misfit  in  our  result  can  be  attributed  to  poorly  constrained, 
shallow  Earth  structure.  The  use  of  shorter  periods  (25-125s)  introduces  a  potential  for 
contamination  from  larger  errors  due  to  unresolved  wave  speed  structures.  At  shorter  periods, 
Earth  structure  will  have  a  greater  influence  on  the  propagating  wavefield  and  errors  in  the  Earth 
model  may  be  mapped  into  the  source  solution.  There  is  a  small,  but  significant,  increase  in 
standard  error  on  the  trend  and  plunge  of  the  P-axis  (or  T-axis)  as  well  as  a  decrease  in  the 
average  cross-correlation  value  from  60-125s  to  25-125s  (Figure  7).  We  do  not  feel  the 
degradation  in  solution  stability,  as  seen  by  the  increase  in  trend  and  plunge  standard  error,  and 
wavefonn  metrics,  as  seen  by  the  decrease  in  the  cross-correlation  value,  are  justification  for  the 
removal  of  shorter  period  signal.  Alternatively,  the  addition  of  shorter  period  signals,  25 -125  s, 
allows  for  much  more  seismic  data  to  be  incorporated  into  the  inversion  indicating  a  robust 
predictive  capability  of  the  3D  wave  speed  model  over  a  wide  period  range.  An  azimuthal  gap  in 
seismic  stations  seen  in  Kazakhstan  and  Russia  also  contributes  to  the  error  of  the  trend  of  some 
solutions,  but  for  most  events  the  overall  azimuthal  coverage  is  excellent  and  the  trend  is  well 
constrained  for  the  entire  data  set  as  a  whole  when  using  the  3D  moment  tensor  inversion. 
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Figure  7.  Standard  Error  for  the  P-axis  of  the  Moment  Tensor 

Event  depth  is  an  essential  variable  when  undertaking  the  adjoint  tomographic  inversion. 
Our  results  for  depth  when  compared  to  the  original  solution  depth  are  shown  in  Figure  8.  Minor 
changes  in  depth,  especially  for  shallow  events,  can  severely  modify  the  resulting  waveforms, 
and  errors  in  the  simulated  depth  can  be  mapped  directly  into  errors  in  the  resulting  3D,  adjoint 
wave  speed  model.  If  errors  in  depth  are  allowed  to  be  mapped  into  the  adjoint  model  they  can 
substantially  affect  the  absolute  wave  speeds,  amplitude  predictions,  attenuation  measurements, 
and  result  in  larger  standard  errors  of  the  final  model.  In  Figure  8  we  compare  our  inversion  for 
depth  with  the  original  depth.  Our  results  compare  quite  well  when  using  a  3D  wave  speed 
model  at  all  period  ranges.  As  with  the  P-axes,  we  performed  a  bootstrap  analysis  to  determine 
the  standard  error  for  the  event  depth.  For  all  events,  the  standard  error  is  smaller  than  the 
symbol  size  in  Figure  8.  Results  at  periods  between  25-125  s,  higher  frequency,  are  more 
consistent  with  the  original  depth  than  those  at  60  -  125s,  low  frequency.  There  is  some 
discrepancy  at  depth  around  100  km  for  the  longer  period  solutions.  This  variability  might  be 
due  to  how  the  initial  solutions  are  constructed  and/or  a  more  limited  depth  resolution  at  longer 
periods.  Variability  is  also  apparent  at  the  shallowest  depths.  The  initial  solution  has  a  depth 
minimum  at  1 5  km  whereas  our  solutions  are  allowed  to  vary  throughout  the  crust.  As  expected, 
many  of  the  shallowest  events  that  were  originally  set  at  1 5  km  have  been  revised  to  shallower 
depths.  Further,  the  use  of  shorter  period  arrivals  produces  more  consistent  depths  than  the 
longer  periods  solutions. 
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Figure  8.  Comparison  Between  Original  and  Update  MT  Solution  Depth 


3.2.  Adjoint  Tomographic  Inversion 

Improvements  in  the  wave-speed  model  for  a  sample  set  of  iterations  is  shown  in  Figure  9  as 
changes  in  the  traveltime  and  in  Figure  10  as  absolute  shear  wave  speeds,  Vsv  and  Vsh.  The 
incorporation  of  transverse  isotropy,  also  known  as  polar  anisotropy,  Vsv  and  Vsh,  is  essential 
for  solution  stability  in  the  mantle.  Figure  9  demonstrates  the  improvement  in  the  wave  speed 
model  manifested  as  misfit  between  data  and  synthetics  seismograms.  The  gray  histograms  are  a 
previous  inversion  and  the  colored  histograms  are  the  improved  model.  Identification  of  body 
wave  and  surface  waves  are  defined  by  the  absolute  wave  speed  of  the  arrival  at  a  station. 
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Figure  9.  Improvements  to  the  Traveltime,  dT,  Metric  at  Select  Model  Iterations 
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A  majority  of  the  improvement  in  the  model  is  on  the  vertical  and  radial  components  as 
measured  by  traveltime,  indicating  an  larger  improvement  to  the  vertical  shear  wave  speed  in  the 
mantle  over  the  horizontal  shear  wave  speed.  Improvements  are  equally  divided  between  body 
and  surface  waves  demonstrating  a  well-balanced  inversion  procedure  that  does  not  emphasize 
one  type  of  arrival  over  another.  This  is  primarily  a  result  the  adjoint  source  creation  that 
normalizes  the  amplitude  of  arrivals  and  relies  on  the  relative  misfit  of  an  arrival  rather  than  the 
absolute  amplitude.  Moreover,  the  inversion  is  dominated  by  misfits  in  traveltime  over 
amplitude  misfit.  This  combination  of  the  body  and  surface  waves  should  produce  a  more  robust 


iteration  model  8 


Figure  10.  Comparison  of  the  Middle  East  Wave-Speed  Model  at  Different  Iterations 
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model  that  has  a  greater  predictive  power  over  a  3D  wave  speed  model  derived  from  a  single 
data  source.  Along  with  the  improvement  of  traveltime  metrics,  amplitude  misfits  are  slightly 
improved.  This  minor  improvement  indicates  that  either  the  bending  of  wavefronts  at  the  current 
minimum  period  are  not  important  or  the  long-wavelength,  anelastic  attenuation  structure  is 
reasonably  characterized  in  the  model.  The  progression  to  shorter  periods  will  begin  to 
incorporate  both  wavefront  bending  and  anelastic  attenuation  in  a  more  dramatic  fashion.  As  the 
iterations  continue  to  localize  wave  speed  structures,  the  possibility  of  severe  wavefront  bending 
becomes  a  greater  possibility.  Moreover,  as  even-more  shorter  period  arrivals  are  incorporated 
into  the  inversion,  anelastic  attenuation  and  variations  within  will  more  dramatically  affect  the 
amplitudes  of  these  arrivals,  necessitating  the  need  to  improve  the  anelastic  attenuation  model. 

As  the  number  of  iterations  increases,  improvements  are  applied  to  the  model  and  shorter 
period  data  are  included  in  the  inversion.  With  these  improvements  and  shorter  period  data 
higher  resolution  structures  begin  to  appear  in  the  wave-speed  model.  Two  examples  of  new 
higher  resolution  structures  are  in  the  Red  Sea  in  Vsv  and  in  the  South  Caspian  Sea  in  Vsh.  It  is 
interesting  to  note  the  wave-speed  model  prefers  reductions  in  wave  speed  in  Vsv  and  increases 
in  wave  speed  as  Vsh;  a  feature  that  would  not  be  possible  without  transverse  isotropy  in  the 
mantle  and  the  fitting  of  both  Love  and  Rayleigh  waves. 

The  incorporation  of  anisotropy  was  essential  for  stability  of  the  inversion  for  wave 
speeds  in  the  upper  mantle.  The  simplest  parameterization  of  transverse  isotropy  with  a  vertical 
symmetry  axis  allowed  for  inversion  stability  and  imaging  of  the  large  wavelength  structure  of 
the  ocean  and  continent.  Inclusion  of  shorter  period  arrivals  will  begin  to  show  the  limitations 
of  this  anisotropic  parameterization.  Work  by  [9]  indicates  anisotropy  beneath  the  ocean  that  is 
well  characterized  by  global  flow  models  except  in  regions  of  upwelling  and  downwelling.  In 
contrast,  the  continents  require  an  understanding  of  the  geologic  strain  history  to  accurately 
define  the  total  amount  of  strain  and  thus  the  anisotropy.  In  essence,  the  oceanic  mantle  can  be 
well-characterized  by  horizontal  flow  and  a  transverse  isotropy  with  a  vertical  symmetry  axis, 
but  the  continents  and  regions  with  vertical  mantle  flow  require  a  more  complex  treatment  of 
anisotropy  appropriate  for  the  direction  of  mantle  flow  or,  in  the  case  of  the  Middle  East,  the 
preserved  strain  history  of  the  Asian,  Arabian,  and  Indian  tectonic  plates.  Without  knowing  the 
true  lateral  changes  in  anisotropy  and  including  large  continental  provinces  within  the  model  it  is 
most  appropriate  to  parameterize  the  mantle  with  a  more  general  fonn  of  anisotropy  to  accurately 
resolve  the  wave  speed  structures. 

A  set  of  comparisons,  waveforms  and  cross-sections,  between  the  current  iteration  and 
the  starting  model  (S2.9EA)  are  contained  in  Figure  11.  Wavefonn  fits  from  a  path  from  the 
southern  Zagros  to  station  KBK  demonstrates  that  updates  to  the  model  are  improving  the  fit 
between  data  and  synthetics.  As  mentioned  before,  this  event  is  not  included  in  the  adjoint 
tomography,  but  is  used  to  provide  a  consistency  check  on  the  model  improvement  and 
progression.  After  about  a  third  of  the  total  iterations,  a  lower  wave  speed  Vsv  lithosphere  is 
required  in  the  middle  of  this  path,  to  the  southwest  of  Tibet’s  western  syntaxis.  This  lower 
wave-speed  region  appears  in  the  cross  sections,  but  is  also  apparent  in  the  synthetics.  The 
surface  wave  in  the  starting  model  arrives  earlier  than  the  data,  but  arrives  “on-time”  using  the 
current  iteration  model. 

Updates  to  the  model  also  include  density  and  compressional  components.  The  initial  3D 
compressional  wave-speed  model  was  scaled  in  the  initial  model  from  the  3D  shear  wave  speed 
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by  a  constant  factor.  Shear  wave  speeds  in  S2.9EA  were  derived  from  a  large  number  of  Love, 
Rayleigh,  and  shear  body  wave  arrivals.  A  comparison  between  the  compressional  adjoint 
results  and  the  initial  model  can  be  difficult  as  the  ratio  between  shear  and  compressional  wave 
speeds  can  vary  due  to  composition.  All  adjoint  iterations  for  compressional  wave  speeds  show  a 
much  slower  upper  mantle  beneath  much  of  the  Middle  East  than  indicated  by  S2.9EA.  Updates 
to  the  model  at  long  wavelengths  for  these  initial  iterations  are  predominately  controlled  by 
Rayleigh  waves  (shear  and  compressional  surface  wave)  as  the  initial  S2.9EA  model  fit  the  Love 
waves  (shear  only  surface  waves)  better  (Figure  1 1). 


Figure  11.  Demonstration  of  the  Current  Iteration  Model  Through  a  Single  Station  -  Event 

Pair 
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Figure  12.  Absolute  Vsv  Speeds,  Initial  vs  Final  Model,  Presented  at  Selected  Depths 
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Figure  12  shows  a  comparison  between  the  Middle  East  3D  wave  speed  model  from  the 
initial  model  to  the  final  model.  Shown  is  the  vertical  shear  wave  speed,  which  is  sensitive  to 
misfits  from  Rayleigh  waves  and  to  a  lesser  extent  body  waves.  The  overall  pattern  of  the 
anomalies  does  not  change  from  the  initial  to  the  final  model,  but  the  absolute  magnitude  and  the 
shape  of  the  anomaly  is  significantly  changed.  As  the  iterations  progress,  anomalies  spatially 
shrink,  but  increase  in  overall  magnitude,  thus  localizing  the  structures.  Examples  of  this 
structure  localization  can  be  identified  at  50  km  in  western  Tibet  and  south  of  the  Black  Sea  at 
100  km.  Further,  artifacts  from  the  initial  3D  wave  speed  model  due  to  either  poor  sampling  in  a 
region  or  ambiguous  structures  not  resolvable  using  traditional  inversion  techniques  are 
diminished  in  amplitude.  Examples  of  artifact  removal  are  at  50  and  100  km  to  the  west  of  the 
Persian  Gulf. 


4.  CONCLUSIONS 

Adjoint  iterations  to  update  the  wave-speed  model  of  the  Middle  East  have  shown  dramatic 
improvements  based  on  data  /  synthetic  comparisons  and  expected  behavior  when  introducing  a 
finite-frequency  tomographic  approach.  Before  iterations  began,  a  seismic  wavefonn  database 
with  -200  reinterpreted  sources  in  the  Middle  East  was  compiled  to  avoid  the  mapping  of  source 
errors  into  wave-speed  structure  [5],  All  sources  in  the  database  agree  with  previously  published 
solutions  and  between  the  methodologies  used:  CAP  [9],  teleseismic  [10],  and  a  MT  inversion 
[2].  Initial  comparisons  between  data  and  synthetics  computed  from  the  initial  3D  model  S2.9EA 
[4]  at  25  seconds  and  longer  show  reasonable  agreement  and  a  general  improvement  over  the  ID 
PREM  model  [3], 

To  further  improve  the  fit  between  data  and  synthetics,  the  adjoint  tomographic  inversion 
methodology  was  implemented  to  iteratively  update  the  3D  wave-speed  model  [1].  The  wave- 
speed  model  was  parameterized  to  incorporate  transverse  isotropy  in  the  mantle  to  accommodate 
the  simultaneous  fitting  of  Rayleigh  and  Love  waves  and  stabilize  the  inversion  in  the  upper 
mantle.  As  the  model  progressed  through  a  series  of  iterations,  regions  of  wave-speed  anomalies 
increased  in  strength  from  the  initial  to  the  current  iteration  as  physics  of  wave  propagation  were 
incorporated  into  the  inversion;  see,  e.g.,  the  Red  Sea  and  the  southern  Caspian  Sea,  Figures  10 
and  12.  Compressional  wave-speed  estimates  also  improved  from  increased  use  of  Rayleigh 
waves.  The  final  model  fits  body  waves  down  to  a  minimum  period  of  14  seconds  and  surface 
waves  to  a  minimum  period  of  20  seconds,  and  the  wave-speed  model  is  able  to  resolve  wave- 
speed  features  with  scale  lengths  of  60  km.  Our  recommendations  for  a  robustly  constrained  3D 
wave-speed  model  in  the  Middle  East  include  using  quality  controlled  seismic  waveform  data 
from  well-constrained  sources  and  use  of  a  physics  based  tomographic  technique  to  improve 
upon  existing  3D  models.  Further  recommendations  to  improve  the  wave-speed  model  include 
adding: 

•  anisotropic  effects  where  necessary,  e.g.  the  mantle,  and 

•  attenuation  to  improve  fits  to  amplitudes,  and 

•  shorter  period  signals  to  further  improve  the  waveform  predictive  capabilities  of  the 
wave-speed  model  of  the  Middle  East. 
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